A Nonlinear ODE System for the Unsteady Hydrodynamic Force – A New Approach

We propose a reduced-order model for the instantaneous hydrodynamic force on a cylinder. The model consists of a system of two ordinary differential equations (ODEs), which can be integrated in time to yield very accurate histories of the resultant force and its direction. In contrast to several existing models, the proposed model considers the actual (total) hydrodynamic force rather than its perpendicular or parallel projection (the lift and drag), and captures the complete force rather than the oscillatory part only. We study and provide descriptions of the relationship between the model parameters, evaluated utilizing results from numerical simulations, and the Reynolds number so that the model can be used at any arbitrary value within the considered range of 100 to 500 to provide accurate representation of the force without the need to perform timeconsuming simulations and solving the partial differential equations (PDEs) governing the flow field. Keywords—reduced-order model, wake oscillator, nonlinear, ODE system I. PROBLEM DESCRIPTION W HEN a uniform flow is interrupted by an infinitelength cylinder, whose axis is perpendicular to the flow, there is a threshold of undisturbed velocity (for a certain fluid viscosity and cylinder diameter) over which a Hopf bifurcation occurs and the steady symmetric wake becomes periodic due to the existence of alternating vortices being shed in the wake at a nondimensional frequency, using the cylinder diameter and undisturbed velocity, that is known as the Strouhal number (St). Experimental and numerical studies showed that this threshold corresponds to a nondimensional velocity, using the cylinder diameter and kinematic fluid viscosity, that is known as the Reynolds number (Re) near 50 [1]–[3]. Before the bifurcation, the resultant hydrodynamic force on the cylinder, due to surface pressure and shear stresses, is parallel to the undisturbed flow. When the threshold Reynolds number is exceeded, this force becomes alternating in both the amplitude and direction, but its average has a non-zero value in the parallel direction. It is a common practice in hydrodynamics and aerodynamics to decompose this force into two orthogonal projections and represent them as a ‘nondimensional’ lift coefficient CL and ‘nondimensional’ drag coefficient CD by scaling these projections with a reference force (per unit length) that is the product of the fluid density, square of the undisturbed velocity, and the cylinder diameter. In this study, we consider the total force directly but we still use a nondimensional representation, which we denote as the total-force coefficient (CT ). We denote the angle of the total force relative to the undisturbed-flow direction (being positive in the clockwise direction) by β as illustrated in Fig. 1. Typical e-mail: [email protected]. histories of CT and β (at Re = 300) are presented in Fig. 2. In these figures (and all subsequent ones), the time is normalized using the cylinder diameter and the undisturbed velocity. These histories require performing numerical simulation and solving the incompressible Navier-Stokes equations (continuity and momentum conservation laws), which govern the unsteady velocity and pressure fields of the constant-density, constanttemperature fluid, which is followed by a post-processing step to compute the surface pressure and shear stress at the cylinder surface and then integrating them to find the resultant force. Whereas well-established techniques and numerical schemes have been developed for this purpose [4]–[6], solving this system of nonlinear, coupled partial differential equations with high accuracy requires a lot computational resources and time. It is therefore desirable to circumvent these simulations by replacing this PDE system by a reduced ODE one that still describes accurately CT and β but with much less computational demands and without the need of post-processing of the velocity and pressure fields. This is the main objective of our study. Whereas details about the flow field can not be obtained through the reduced-order model, the force on the cylinder is of primary importance when we are concerned with the design of a cylinder (or a pipe) and the implied fatigue problem. A full cycle of β corresponds to a full cycle of shedding where two contra-rotating vortices are shed in the wake. We marked four equally-spaced instants of time over approximately one cycle of β starting from a point where β is maximum. The corresponding distributions of the surface pressure (nondimensional) at these instants are presented in Fig. 3, and we superimposed an arrow indicating CT and its direction. We should mention here that we only show the distributions of the surface pressure and not the surface shear because the former is very small compared to the latter (except at very low Re before shedding takes place, which is not of interest here) as also indicated in Ref. [7]. However, in all the results of CT and β we obtained through solving the NavierStokes equations and present in this study, we accounted for the contribution of both the pressure and shear. The average (over time) surface pressure distribution is symmetric and the average CT is not zero but the average β is zero as indicated in Fig. 4. For the same case of Re=300, we present in Fig 5 typical spectra of CT and β. We normalize the frequencies by the Strouhal number. As mentioned earlier, the fundamental frequency of β is the shedding frequency. The fundamental frequency of CT is twice the fundamental frequency of β, thus is twice the Strouhal number. The total-force coefficient can be represented as a superposition of an average value (denoted as a0T ) and a harmonic function with a frequency at 2St and World Academy of Science, Engineering and Technology International Journal of Mathematical, Computational, Physical, Electrical and Computer Engineering Vol:4, No:3, 2010 413 International Scholarly and Scientific Research & Innovation 4(3) 2010 scholar.waset.org/1999.7/4282 In te rn at io na l S ci en ce I nd ex , M at he m at ic al a nd C om pu ta tio na l S ci en ce s V ol :4 , N o: 3, 2 01 0 w as et .o rg /P ub lic at io n/ 42 82 amplitude a2T , which is the magnitude of the CT harmonic at a frequency equal to 2St. All odd harmonics and other even harmonics of CT are negligible as indicated in Fig 5. On the other hand, all even harmonics of β (at even multiples of St) are negligible. The magnitude of the third harmonic a3β is one order of magnitude smaller than the fundamental harmonic (whose magnitude is a1β), and the higher odd harmonics are negligible. These features of CT and β are common for the Re range we consider here (from 100 to 500), and we will utilize them when proposing our new ODE system that models CT and β. II. MODELING OVERVIEW AND OUR OBJECTIVES Birkhoff and Zarantonello [8] observed that the wake of a stationary cylinder, with its continuous swinging, can be modeled by an oscillator and alluded to a linear one. Seven years later, Bishop and Hassan [9] conducted a thorough experimental study of a cylinder forced to oscillate perpendicular to a uniform flow. When the forcing frequency was close to the Strouhal number, the wake responded at the forcing frequency and not at the Strouhal number, and we speak of a wake synchronized with the motion of the cylinder. Based on this observation, Bishop and Hassan proposed that the periodic wake can be modeled by a simple oscillator. They did not specify a particular oscillator (which they called ’wake’ or ’fluid’ oscillator), but indicated that it is nonlinear and selfexcited. Since then, several models have been proposed for the wake of a rigid cylinder [10]–[23]. Each model is a nonlinear ODE that, when integrated in time, yields the history of the lift coefficient. Two candidate models have been commonly used so far for the lift force on a cylinder in a uniform flow: the Rayleigh oscillator [10], [11], [13], [14], [19] and the van der Pol oscillator [20]–[23]. The van der Pol oscillator for the lift coefficient (CL) has the form C̈L + ω 2 CL + μ ĊL + c C 2 L ĊL = 0 (1) and the Rayleigh oscillator has the form C̈L + ω 2 CL + μ ĊL + c Ċ 3 L = 0 (2) Both are self-excited self-limiting nonlinear oscillators with a cubic nonlinearity. Krenk and Nielsen [18] proposed a combination of these two oscillators, and Griffin et al. [24] and Skop and Griffin [25] even added a Duffing term (C L). Landl [15] proposed a variant of the van der Pol oscillator with an additional quintic term (C L ĊL) as an additional hydrodynamic damping, but this increased the model parameters to be determined and its complexity. It should be noted that only odd nonlinearities can be considered in the wake oscillator for CL because its spectrum mainly consists of a dominant component at St and harmonics at odd multiples of St. It was Nayfeh et al. [20] who justified their choice of the van der Pol oscillator, rather than the Rayleigh one, to describe the lift coefficient exerted on a stationary cylinder by the wake. They reinforced that the van der Pol oscillator produces a phase of the third harmonic, at 3St, relative to the fundamental one, at St, that is equal to 90 as compared to 270 in the case of the Rayleigh oscillator. They found that corresponding phase obtained from numerical simulations (solving the Reynolds-averaged Navier-Stokes equations at different Re) was approximately 90. Qin [22] and Marzouk et al. [23] based their choice of the same oscillator on the reasoning of Nayfeh et al. [20], whereas Facchinetti et al. [21] did not give a reason for choosing this oscillator. The wake oscillator of Facchinetti et al. [21] lacks the flexibility of describing the changes in the lift coefficient at different wake regimes (different Reynolds numbers) because the coefficient of the (negative) linear damping and (positive) nonlinear damping are equal, thus the model yields a unique limit cycle irrespective of the Reynolds number, which is a serious violation of the physics of the problem. The model also assumes a constant Strouhal number of 0.2, which also should be allowed to vary with the Reynolds number as found in experimental and numerical studies (see Appendix A). We note that the model variable in the work of Facchinetti et al. [21] is the CL scaled to a reference constant and not the CL itself. However, this does not affect the model dynamics. Marzouk et al. [23] added a Duffing term to the van der Pol oscillator. Although Griffin et al. [24] and Skop and Griffin [25] have previously included a Duffing term in their mixed (van der Pol and Rayleigh) wake oscillator, they did not provide a reason behind including such term, which seems to be just an attempt to generalize the oscillator. In contrast, Marzouk et al. [23] explained the need for this term by the ability of the model to capture exactly the relative phase of the third harmonic of CL, which would be fixed at 90 o without this Duffing term. They also provided expressions relating the model parameters to the CL characteristics (e.g., the frequency and magnitude of the fundamental harmonic). Ogink and Metrikine [26] started from where Facchinetti et al. [21] ended and modified the cubic term in the van der Pol oscillator from C L ĊL to 1/(a + b C 2 L)ĊL, where 0< a <1 and 0< b are two tuning parameters. This modification was proposed so that for increasing values of CL, the influence of the nonlinearity in the modified model decreases. However, they did not discuss how a and b are related to the Reynolds number or the wake characteristics. They set these extra tuning parameters to constant values of a=1/2 and b=1. Some studies [16], [20], [22], [23], [27] also considered a reduced-order model for the drag. Currie and Turnbull [16] used a Rayleigh oscillator to describe the oscillatory part only of the drag coefficient. Assuming a harmonic oscillatory drag, a relationship between the parameters of the linear and nonlinear damping was established such that the amplitude of the oscillatory drag coefficient becomes 0.2, and the model does not allow any variation of this value with the Reynolds number. This is a remarkable drawback in the model, which also specifies the individual values of these parameters based on matching the model results to an arbitrary data set in a trial-and-error fashion. Kim and Perkins [27] considered an elastic cable suspension, and introduced a drag model in which the oscillatory part of the drag coefficient was represented as the product of a time-dependent function and a spatial function whose argument is the longitudinal coordinate (divided by the cable World Academy of Science, Engineering and Technology International Journal of Mathematical, Computational, Physical, Electrical and Computer Engineering Vol:4, No:3, 2010 414 International Scholarly and Scientific Research & Innovation 4(3) 2010 scholar.waset.org/1999.7/4282 In te rn at io na l S ci en ce I nd ex , M at he m at ic al a nd C om pu ta tio na l S ci en ce s V ol :4 , N o: 3, 2 01 0 w as et .o rg /P ub lic at io n/ 42 82 diameter) along the cable in the static condition. The timedependent function was modeled by a van der Pol oscillator whose natural frequency is at twice the Strouhal number. A similar oscillator was used for the lift coefficient but the natural frequency is at the Strouhal number. Kim and Perkins used seven quadratic coupling terms in the two oscillators with seven parameters that need to be identified. These parameters were grouped in four functions, whose values were simply selected manually to fit some experimental data, which also was followed in determining other model parameters for the uncoupled part of the oscillators. Therefore, their coupled lift and drag models have many parameters to identify; there is no robust method of determining these parameters, and the average drag is not reproduced. Nayfeh et al. [20], Qin [22], and Marzouk et al. [23] assumed the oscillatory part of the drag coefficient to be a harmonic function at twice the Strouhal number and modeled it by an algebraic function of the independently-modeled lift coefficient. Nayfeh et al. [20] used a single-term quadratic function (proportional to CL ĊL). This freezes the phase between the drag and lift to 270. Qin [22] used another single-term quadratic function (proportional to the oscillatory part of C L) and added a linear coupling term (proportional to CL). However, our simulations show that such linear coupling (which was proposed only by Qin) is not necessary. Marzouk et al. [23] extended the drag model of Nayfeh et al. [20] by combining their coupling term with the one used by Qin [22], so that the resulting two-parameter algebraic model can reproduce any value of the relative drag phase. The coefficients of these terms were analytically related to the relative drag phase. Marzouk and Nayfeh [28] modeled the total-force coefficient and its angle, and considered a mixed van der Pol and Duffing nonlinear terms for the angle oscillator, but used an algebraic quadratic equation to relate the oscillatory part of the total-force coefficient to the evolving angle. Thus, their differential/algebraic wake model cannot resolve the average value of the total-force coefficient. In addition, they considered a single Reynolds number. The present study addresses these limitations. The quality of a reduced-order model for the force on a cylinder due to the periodic distributions of the pressure and shear stresses on its surface as a result of the nearwake shedding process is gauged by the capability of the model to capture qualitatively and quantitatively the main characteristics of this force, which inevitably depend on the Reynolds number and not universal constants. Therefore, we propose in this study a new wake model, which consists of two coupled ODEs for the nondimensional (or coefficient CT ) resultant force of the instantaneous surface stresses and its angle. The model has three advantages: First, it directly models the actual hydrodynamic force rather than modeling the lift and drag, which are two ‘convenient’ components of this actual force. It is still possible to retrieve the lift and drag in a simple post-processing step of the resolved total force, where CL = CT sin(β) and CD = CT cos(β). Second, we provide closed-form expressions relating the parameters of the model to the main flow characteristics, which depend on the Reynolds number. Therefore, our ODE system implicitly accounts for the variations of the wake regime with the Reynolds number and does not freeze any of these characteristics (such as the Strouhal number or amplitude of the lift coefficient). The model parameters can be obtained in a schematic way which is faster and more accurate than manually matching these parameters to match a target set of data. Third, our proposed ODE system captures directly the average component of the total force (thus the average drag) and does not ignore it. This, in addition to directly modeling the actual hydrodynamic force, makes the model better than the existing ones in terms of being able to mathematically describe the physical aspects of the problem with simple equations which are much easier to solve than the full Navier-Stoked equations. We provide continuous functions of the model parameters, within the considered range of Reynolds number, so that one can use the proposed model at any arbitrary value of Reynolds number, which can then be coupled with a structural oscillator to investigate the two-degree-of-freedom fluid-structure interaction between the elastically-mounted cylinder and the flow. The modeling concept we followed here is not strictly limited to a circular cylinder. Rather, it can be extended to other geometries that exhibit a similar shedding phenomenon [29]. III. PROPOSED ODE SYSTEM As indicated in the simulation results in Fig. 5, the spectrum of β consists of a dominant component at St (with magnitude a1β) and a third harmonic at 3St (with magnitude a3β). Such behavior of β suggests a self-excited oscillator such as the van der Pol or the Rayleigh oscillators. Looking at the phase by which of the component at 3St leads the one at St (which we denote by ψβ3) over the range of Re under consideration, we found it to lie between 148.8 and 154.5. Therefore, the van der Pol oscillator (which implies a modeled ψβ3=90 ) is preferred to the Rayleigh oscillator (which implies a modeled ψβ3=270 ). However, we still need to adjust the van der Pol oscillator by either combining it with a Rayleigh oscillator or a Duffing one, so that the weights of the two cubic-nonlinearity terms in the combined oscillators can be adjusted so that the modeled ψβ3 is precisely equal to a desired value that is known to occur at a certain Re. Whereas for the case of modeling the force on a stationary cylinder, there is no obvious preference of which oscillator to combine with the van der Pol oscillator, we examine the case of a cylinder undergoing harmonic oscillations to make this choice. Experimental [9], [30], [31] and numerical [32]–[34] studies at low and high Reynolds numbers show that the frequency-response curves of the CL, CD, and surface pressure are asymmetric about the St. When the oscillation amplitude is high enough, hysteresis and jumps occur, which coincide with a transition between two different wake modes (one having stronger forces than the other). The forced van der Pol, Rayleigh, or a combination exhibits symmetric frequency-response curves about the linear frequency, and two symmetrically-located hysteresis locations would be observed in that case [35], which does not match the physics of the problem. On the other hand, adding the Duffing term to the van der Pol oscillator can produce biased World Academy of Science, Engineering and Technology International Journal of Mathematical, Computational, Physical, Electrical and Computer Engineering Vol:4, No:3, 2010 415 International Scholarly and Scientific Research & Innovation 4(3) 2010 scholar.waset.org/1999.7/4282 In te rn at io na l S ci en ce I nd ex , M at he m at ic al a nd C om pu ta tio na l S ci en ce s V ol :4 , N o: 3, 2 01 0 w as et .o rg /P ub lic at io n/ 42 82 frequency-response curves with single location of hysteresis [36] because the backbone curve is no longer vertical (as in the case of forced combination of the van der Pol or Rayleigh oscillators). Therefore, we propose a combined (van der Pol and Duffing) oscillator for β. It is desired to use a linear oscillator for the amplitude of the total-force coefficient CT because of its simplicity and to reduce the number of overall-model parameters. However, proper forcing is needed (unlike the free β oscillator) in order to obtain a stable nontrivialCT response. A functional form of β is a good candidate for this purpose because physically it is the complementary feature needed with CT to fully describe the evolution of the total force coefficient. Utilizing the twoto-one frequency relationship between CT and β, we set the linear frequency of the CT oscillator to be twice the frequency of the β oscillator. We seek a forcing of the CT oscillator that satisfies three conditions. First, it is at a frequency of 2St (which is the correct frequency of the CT ). Second, it has an average value, which provides the average CT . Third, it has at least two terms, so that the resultant phase by which CT leads β is not a universal constant and can vary with the Reynolds number. Under these conditions, the simplest form of the forcing is a combination of β and β β̇. We note that β could be replaced by β̇, but we will use the latter, mainly for convenience as β is a more ‘visible’ quantity than β̇. Based on this discussion and justifications, we propose the following system to describe the time-dependent total-force coefficient (CT ) and its angle, in radians, (β): β̈ + ω β + μ1 β̇ + c1 β 2 β̇ + c2 β 3 = 0 (3a) C̈T + 4ω 2 CT + μ2 ĊT = q1 β 2 + q2 β̇ β (3b) The proposed system in Eq. (3) has seven model parameters; namely the linear frequency ω, the linear ‘destabilizing’ damping coefficient μ1 <0, the cubic ‘stabilizing’ damping coefficient c1 >0, the Duffing-term coefficient c2, the linear ‘stabilizing’ damping coefficient μ2 >0, the quadratic static-static coupling coefficient q1, and the quadratic statickinematic coupling coefficient q2. In the remaining part of this section, we analyze the model and derive expressions for these parameters, which analytically relate them to seven target flow characteristics, which we assume they are known (e.g., from processing of a simulated or measured flow field or from interpolating the parameters from pre-determined values at different Reynolds numbers). These flow characteristics are the Strouhal number (St); the harmonics magnitudes a1β , a3β , a2T ; and the relative phases ψβ3 and ψT2. We use the method of harmonic balance and find first the algebraic equations for the four parameters in Eq. (3a) and then for the three remaining parameters in Eq. (3b). In the following analysis, we do not make any assumptions about the relative magnitudes of these parameters. A general Fourier series for the periodic β, with the dominant component and its small superharmonic component can be written as β (t) = a1β cos(2π St t+ σβ) + a3β cos(6π St t+ 3σβ + ψβ3) (4) Because β has a self-excited limit cycle and no external forcing appears in Eq. (3a), σβ1 is arbitrary. So, we can set σβ1=0. With this, we re-write Eq. (4) as β (t) = a1β cos(2π St t) +a3β cos(6π St t) cos(ψβ3) −a3β sin(6π St t) sin(ψβ3) (5) Substituting Eq. (5) into Eq. (3a), and then collecting the coefficients of cos(2π St t), sin(2π St t), cos(6π St t), and sin(6π St t) results in the following system of algebraic equations: χ1 ω 2 + χ3 c1 + χ4 c2 = χ5 (6a) δ2 μ1 + δ3 c1 + δ4 c2 = 0 (6b) 1 ω 2 + 2 μ1 + 3 c1 + 4 c2 = 5 (6c) γ1 ω 2 + γ2 μ1 + γ3 c1 + γ4 c2 = γ5 (6d)